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1. introduction 


Recent aovances in the computation of transonic potential flows 
around airfoils and around complex aircraft have provided powerful 
tools for the aerodynamic design of an aircraft (Caughey and 
Jameson, 1979; Mercer and Murman, 1980). These computational codes 
are used extensively in practical aircraft designing to the extent 
that the entire aerodynamic configuration can be determined by using 
computational techniques (Jupp, 1980). These advances can also be 
readily applied to the flow around a rotor in turbomachinery. 

The applicability of the steady potential flow -assumption to 
flow around a rotor has been discussed in the past (Jou, 1979; 
Dulikravich, 1979). The main conclusion is that in a rotating 
fixed-blade frame of reference, a reduced potential exists such that 
its gradient is the perturbation velocity. Both Dulikravich (1979) 
and Adamczyk (1980) have applied the finite volume technique to 
flows around cascades, and Dulikravich has extended the computations 
to three-dimensional flows. The features of these two works are 
listed in the following table. 



Dulikravich 

Adamczyk 

Mesh 

0-type: Conformal mapping 
and shearing 

0-type: Generated by 
an electrostatic 
problem 

Full Potential 

Yes 

Yes 

Conservative 

Yes 

Yes 

Artificial 

Dissipation 

Artificial Viscosity 

Artificial Density 

Solution 

Algorithm 

SLUR 

SLOH 

Dimensions 

Two and Three 

T wo 
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There are a few questions remaining to be answered. Due to geo- 
metrical constraints, the 0-type mesh cannot sufficiently track a 
vortex wake far downstream. Whether the induced flow is accurately 
computed in the three-dimensional case remains to be verified. A 
C-type mesh seems to be more desirable in the three-dimensional com- 
putations. The different methods of implementing the artificial 
dissipation in the artificial viscosity formulation and the artifi- 
cial density formulation (Hafez et al., 1978) need a fair comparison 
in their merits. This can be achieved only by using a code with 
exactly the same algorithm except for the dissipation terms. The 
last and one of the most important questions is whether an aigoiithm 
more efficient than the successive line overrelaxation scheme can be 
developed. An efficient algorithm is an important element in using 
the numerical method for design purposes. 

Tlie present work focuses on the following aspects of the finite 
volume calculation; 

(1) Development of a two-dimensional finite volume SLOR code 
using a C-type mesh coordinate system. 

(2) Development of a scheme with both the artificial viscosity 
and artificial density options and comparison of the 
results of the computation. 

(3) Investigation of possible means to accelerate ttie 
convergence of the computation. 

(A) Recommendations for the development of a three-dimensional 
code based on the two-dimensional investigation. 


Dr. John Auaniczyk has provided the original mesh yeneratiori 
scheme PJUHID tor modification in this work. 
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2. Finite Volume Method. Artificial Density Formulation, and 
Nonlinear Iteration 


The finite volume method has been reportei in detail by Jameson 
and Caughey (1977) and does not require an extensive review here. 
However, there is some confusion about (1) whether the method is a 
finite difference or finite element method, (2) how the nonlinear 
term is linearized during the iteration, and (3) how an artificial 
density concept would be implemented. The present section intends 
to clarify these que^'ions. 

The finite volume scheniti as implemented ty Jameson and Caughey 
(1977) consists of the following procedure: 

(a) Mapping of the physical domain into a rectangular computa- 
tional domain (construction of curvilinear coordinates). 

(b) Computing the flux phu^ at the center of the cell formed 
by coordinate lines as shown in Figure 1. The flux is 
computed using the local bilinear variation of the potential 
♦ (x^) and the transformation x^^(Xj) . 

(c) Adding a recoupling term which effectively relocates the 

point of evaluated flux phU a distance 6 away from the 
center as shown in Figure 1. Here, 6 < 1/2 . (In the 

present version, the variation of h has not been included.) 



Figure 1 . Local Coordinate Trantformadon 
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(d) Forming the flux balance for the “secondary" cell 
represented by the dotted line in Figure 1. 

(e) Adding the artificial time-dependent term. 

(f) Adding the artificial viscosity term. 

To see the relation between the above procedure and the finite 
element method, one should start with the Bateman variation 
principle; i.e., conservation of mass can be recovered if the 
variation of the following integral is carried out: 


I = / p dfi 



(1) 

Y 

where p = p ' 

1 


(2) 

and p = ^1 + 

1 ^ 

• 

(3) 

By transforming the domain of integration 

from the 

physical 

space 

fi(x) to the computational space the 

integral 

becomes 


I = y* ph dfi 

<1 


(4) 

with p = [l + M^il - 

L ^ 3X 


• 

(5) 


The finite element method amounts to approximating 1 by a sum of 
integrals over small discretized elements: 



n=l F.E. 


To evaluate the integral inside the elements, certain trial functions 
for and x^(Xj) must be assumed. If bilinear dependence of 4> 
on Xj is chosen for the elements, these integrals can be evaluated 
explicitly in terms of the values of 4> at the corners of the 
elements. The variation of the integral with respect to these corner 
values gives the discretized system of equations for 4>. It can be 
shown by fairly tedius calculations that if the integral in the 
element is evaluated by expanding the integrand in a Taylor series 
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around the center of the element, the variational procedure leads to 
a discretized problem corresponding to the finite volume method with 
the displacement factor 6 = 1/6 . Aside from the artificial timo- 
dependent term and the artificial viscosity term, the finite volume 
method is indeed an approximation to the finite element method. 
However, the artificial viscosity term requires evaluations of the 
higher order ^.erivatives , etc.) within the element. Hence, 

with the bilinear element, there is no way that the information of 
artificial viscosity can be contained inside the element. Uo to now 
this term has been treated in the finite difference manner in all the 
existing finite element works (e.g., Deconinck and Hirsch, i9BU). It 
is also the root of the problem in using more flexible triangular- 
type finite element calculations. 

The artificial density version of the dissipation in transonic 
computations (Hafez et al., 1978) has gained some popularicy. In an 
arbitrary mesh, the following derivation gives the proper form of 
the artificial density: 


or 


P 


s 




p = p - (jPgAs 


1 

37T 


1 

T7T 


[up^ * 


- (UV) 


9x 3x 
9 Y 

9X 9Y 



(7) 

( 0 ) 


(9) 


( llJ) 


where (C^^, ) art- the covariant components of tlu v( locily 

vector. Letting ttie norm ot the vector (l., , he netiiieo .r. 

( 1 1 ^ I + I C Y 1 ) , we liave 

“s''" - ‘'"’x * ''“y> ■ 

hq 


Ui) 
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This expression avoids taking the square root in the computation. 

The derivatives ond py are evaluated in the upwind direction, 
while the switching function p is evaluated one-half mesh upwind. 

The use of artificial density is usually coupled with an iterative 
procedure different from that used in the artificial viscosity version. 
In the fully implicit iteration procedure of the artificial viscosity 
version originally formulated by Jameson, the nonlinearity of the 
governing equation is linearized in the following manner; 

= 0 (12) 




_n+i n „ n 
P = P + v^p «vc 


(13) 

(U) 


where C is the correction to the velocity potential and is the 
gradient in the velocity space. Hence, the equation for the correc- 
tion C is 

V'fp'^VC + (V^p^*VC) = -V(p^V4>^) . U5) 


Although Jameson does not take full account of the left-hand side, 
he does include the factors [(1 - u /a ), etc.], which appear in 
the quasi-linear version of tiie governing equation. The iteration 
is, although not to the full extent, an implicit method. The itera- 
tion procedure which is often used in conjunction with the artifi- 
cial density version is 

V*(p'*'7<|)'^*^) = 0 . (16) 

The artificial density o’ is evaluated from the last iteration. Iht: 
The procedure is seiuiimpl ici t . This explicit use of p caused iiali.v 
(1979) to call it an "elliptic” method. The u^^e of the tei miiirjluijy 
is unfortunately not justified. The mattiematical nature of tin; t.(^ua- 
tion can only be defined by discussing the equivalent time-dependent 
problem in which the steady-state problem is embedded. 
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To illuminate this, we consider a simpler model equation (small 
disturbance) ; 


*J'vv ” ® 


p = 1 - 


YY 


(17) 


Now, 


^n „n + l ^n + 1.. 

P = p - p^ At 


(18) 


n + 1 


n + 1 ^ 
'Xt 2 


By substituting these relations into the governing equation with 
a semiimplicit iteration procedure, one produces the following time- 
dependent problem: 

At 


(P4>x)x * * T l^t*^X^X = ° 


ll - i:»X ^ 2 '^xt^'^'xx *^YY ^ 2 *^X^^XXt " ° 


At 


(19) 

( 20 ) 


Since all variables are in the (r» ♦ 1) level, the superscripts are 
omitted. There is an additional artificial time-dependent term 
which is added to enhance the stability (Jameson, 197A). The 
complete time-dependent problem with the artificial dissipation can 
be represented by the following equation: 


-o^Xt * ^*^XX * *^YY * ^"^XXX *^^XXt ■ ° 

where a, A, u , and c are functions of the solution itself. It 
should be remembered that u = OlAX) and c = 0(At). The third-order 
terms are of singular perturbation type. The effects of tfie addi- 
tional term <txxt unsteady solution remain to lie studiecj. Jl 

cannot be pronounced simply as an elliptic problem ir» the spaliiil 
domain . 
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3 . Numerical Code 


The nu't.j*' i.cal code developed here is based on the two- 
dimensional finite volume code, FLO 26, developed by Jameson for 
flow around an airfoil. The modifications include; (1) generation 
of a C-type grid for cascades, (2) the boundary conditions on both 
the periodic boundaries and the upstream and downstream boundaries, 
and (3) the artificial density option in the solver. These are 
discussed in the following subsections, 

3.1 Mesh Generation 


The mesh generation routine in the present work is the modifica- 
tion of an 0-type mesh generator by Adamczyk (1930). The method is 
based on an electrostatic problem of a cascade charged with V = 1 
and V = -1 on alternate blades. The equipotential lines and the 
flux lines are used cs the coordinate lines. To aoopt the method 
for generating a C-type mesh to be used in a finite volume calcula- 
tion, the following modifications are made to the original mesh 
generation routine; 

(1) A straight line is attached to the trailing edge of the 
blade. This straight line is an approximation of the 
trajectory of the vortex wake (in the three-dimensional 
case) and is assumed to be at the angle of the bisecting 
line of the trailing edge. The blade and the straight line 
are charged with the same potential in the electrostatic 
problem. 

(2) The resulting (J-type mesh system around the combined 
geometry of blade and wake is truncated near the Ocwn.streain 
end of trie straight line (the wake). Ttit resulting me^i; is 
of the C type as far as the original blade system is 
concerned . 

(3) To ensure ttiat a niesh line will pass tliruuyh the tiaiiing 
edge, a shearing transformation is applied. This shearing 
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w ransformation also adjusts the number of meshes on the 
blade to a multiple of four. The purpose is to perform 
computations on successively refined meshes to increase thu 
rate of convergence. 

(4) The mesh lines intersecting the "woke" are sheared so that 
they are continuous across the wake. 

Other improvements on the mesh generation scheme are also per- 
formed. Most importantly, the input data of blade geometry is in- 
terpolated into a more coarsely spaced approximate blade shape. An 
approximate mesh system is generated. The mesh lines are then 
sheared to the exact blade position. This procedure increases the 
speed of the mesh generation procedure without loss of accuracy in 
defining the blade geometry. The results of the mesh generation are 
shown in Figures 2, 4, and 9. 

3 . 2 Bou.idary Conditions 

On the periodic line, an extended computational domain is 
established by including the next "exterior” mesh lines, as shown by 
the dotted lines in Figure 3. The values of the potential on these 
"exterior lines" are given by the previous iteration from the 
condition of periodicity. These values are used as the Dirichlet 
boundary conditions. 

On the exit boundary, the boundary condition that the flow field 
attains its asymptotic downstream state is imposed. After each 
iTtT-ation, the circulation around the blade f is computed. The far- 
field conditions are computed by the solution of the following 
equations: 


r:os Cij = cos 


('/'/) 

sin Oj c q^ sin 

I' 

H 

('/^) 




* [l ♦ M^U - 

l] 

(24 ) 
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Fine Grid for NACA 0012 
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wherc jno a.^ are the flow angles far upstream and down- 
stream, respectively; and Q 2 are the density and the speed 

far downstream; and H is the pitch distance. Using the farfieid 

velocity vector (q^, 02 )* values of the velocity potential 
on the boundary points are obtained by integration from the neigh- 
boring interior points. These potential values are used as 
Dirichlet conditions for the next iteration. 

As for the inflow boundary, we choose to use the Neumann bound- 
ary condition. Keferrlng to Figure A, the physical velocity com- 
ponents (u, v) at the center of the inflow boundary ceils, marked 

with the cross, are given by the prescribed upstream condition. 
Contravar iant velocities are computed from the given velocity for 
flux balance. 

These boundary conuitions are implemented numerically and they 
have proven to produce a stable Iteration procedure. 

3.3 Artificial Density option 

The artificial density formulation, with an explicit evaluation 
of the artificial density from the previous iteration as discussed 
in Section 2, is coded as an option in the solver subroutine, with 
this option, the merits of one against the other can easily be 
evaluated, further discussions will be given in section n. 
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Numerical Results and Discussion 


Two blade geometries are chosen for developing the numerical 
code. The NACA 0012 configuration is chosen for comparing the re- 
sults to the single airfoil calculations so that the effects of a 
cascade can be seen. Sanz's shock-free compressor blade is chosen 
to validate the code against the design conditions. 

Figures 5 and 6 show the computed pressure distribution on the 
NACA 0012 with a pitch distance of 1.98 C at Mach number 0.75 

= 0*75) with an angle of attack relative to the blades at 
-2 degrees (a = -2°) and at various stagger angles. The effects 
of a cascade are obvious from a comparison of these figures. These 
results have been obtained by using the artificial viscosity 
version. Figure 7 shows the results of computation using the 
artificial density version with the same artificial dissipation 
coefficient as that used in the artificial viscosity version. There 
are ripples ahead of the shock wave and the solution is slowly 
diverging. With the increase of the dissipation coefficient to 1.8 
times that previously used, the solution converges; the results are 
shown in Figure 8. The convergence rate, as well as the quality of 
the shock resolution, is comparable for these two types of 
formulations for artificial dissipation. 

For Sanz's shock-free blade, the preliminary calculations using 
the given data show unreasonable results near the leading edge. A 
close examination of the blade data reveals a lack of data points 
near the crucial leading-edge region. The data around that region 
are filled in by visual interpolations. Figure Iti siiows the n-sultr. 
of computations for the design condition. Tlte results compaie very 
well with the Mach number distribution for the uesign coiuiitiun. 
These results are computed using the artificial density version. 

The computations of the same case using the artificial viscosity 
version with a coefficient of dissipation as large as 3 are shown in 
Figure 11. Notice the bad results near the leading edge. Experi- 
ences in the past have shown that the artificial viscosity version 
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Figure 6. Pressure Distribution (NACA 0012) (2) 
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Figure 8. Calculation Using Artificial Density 
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Figure 10. Pressure Distribution at Design Condition (Sanz Blade) 
(ALPHA Relative to Horizontal Direction) 
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is extremely sensitive to the discontinuity of the mesh distrihu- 
tion, while the artificial density formulation can tolerate some 
mesh irregularity. 

Figures 12, 13 and 1^ show some results of computations of 
Sanz's blade under off-design conditions. The case with 

= 0.711 and a = 32® shows an interesting depression of 
pressure distribution on the upper surface near x/c = 0.2b. The 
question remains whether a weak shock wave is presented there. The 
resolution of the grid system does not allow us to resolve the 
question. For the case where M = 0.78 and a a, 30.01®, the shock 
wave is so strong that the boundary layer may separate and a strong 
V iscoiis-inviscid interaction will modify the pressure aistribution 
in a najor way. 

From these results, the following conclusions can be made: 

U) The analysis code does produce an accurate solution of the 
transonic flow over a cascade as far as comparison with the 
design condition is concerned. 

12) The artificial density formulation is comparable to the 
artificial viscosity formulation in both accuracy and 
convergence . 

(3) The artificial density option can handle the mesh 

discontinuity better than the artificial viscosity version. 
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Figure 1 2. Pressure Distribution at Off • Design Condition ( 1 i 
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5 . Multigrid Method in Two Dimensions 

As mentioned in Section 1, the convergence of the successive line 
overrelaxation scheme is slow. For the computational method to be an 
efficient tool in the design process, the speed of computation in 
achieving certain accuracy must be increased. There are some recent 
progresses in applications of the multigrid method to transonic flow 
problems. Aside from the earlier report of South and Brandt (1977) 
on the small disturbance equation, Arlinger (1980) has reported suc- 
cesses in applications of the method to multielement configurations 
as well as axisymmetric flows. Jameson (1979) in his recent work has 
coupled the ADI scheme with the multigrid method and has shown some 
'-.pectacular reductions in the effort of computation. In the present 
work, we shall explore the possibility and the method of applying the 
multigrid method to the cascade problem. 

It is believed that the successive line overrelaxation scheme is a 
poor smoothing scheme for high-frequency components of the error in 
the multigrid computation. The blame is not on the relaxation process 
itself but on the stretching of the coordinate lines in the airfoil 
problem. It can be shown that if the mesh size Ax in the marching 
direction is small compared to Ay in the line direction, the ampli- 
fication factor for the iteration procedure approaches unity for high- 
frequency components as Ax/Ay approaches zero (South and Brandt, 

1977). Thus, Jameson (1979) was forced to use the ADI scheme as the 
smoothing scheme in his multigrid process. The situation is quite 
different in the cascade problem. The mesh size in the line direction 
is restricted by the gap distance and requires no stretchirtg at all. 
The line overrelaxation method can be an effective smoothing scheme 
in the multigrid process. Furthermore, the number of meshes across 
the gap is usually limited, while it is quite a large number in the 
other direction. In the computational space, the domain of interest 
is a narrow band. The information exchange in the direction of the 
gap is achieved by the line relaxation; but, the information is 
propagated one mesh per iteration in the direction of the passage. 

It is particularly desirable to use the idea of the multigrid in the 
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direction of the passage, l.e., the longer dimension in the computa- 
tional domain. To test the idea, a preliminary effort in that 
direction has been performed in the present work. 

Following Jameson (1979), the following multigrid scheme is form- 
ulated. For the nonlinear equation L^(<^) = 0, the flnegrid approxi- 
mation is obtained by the smoothing scheme. The residue in the 
fine grid is obtained. The correction in the coarser mesh can be 
obtained as 

'•zhiShJ = 

2 h 

where is an interpolation operator from the fine grid cO the 
coarse grid which will filter out the high-frequency component of the 
residue. But, is a nonlinear operator depending on the coarse 
mesh potential fo use the same nonlinear smoothing scheme (SLOR) 

in the same form as that in the fine mesh, we add to both 

sides of the coarse-mcsh correction equation to form 

'"2h^''''h * ‘^2h^ = ‘-2h^'^2h^ 

= 4hl*h> - 

So, after the fine-grid solution is obtained, the right-hand side 
of the above equation is evaluated. The smoothing scheme (SLOR) 
is again used to obtain . To correct the fine-mesh solution, the 
reverse operation is applied; i.e.. 


‘^2h = *2h - 

K1 

L J at 2h mesh 

(27) 

n + 1 .n 

h = *h * 

^2h ^Sh ^ 

(28) 


where n denotes the level of iteration and 1^^ is the interp(jJ at i on 
operator which interpolates the correction from the coarse grid 
to the fine grid. 
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Kor the cascade problem, the following weighted average is used 

O K 

for the interpolation operator at the mesh point (i|j), 



(i) 



[' 






24> 


i] 


(29) 


while the linear interpolation is used for the coarse-to-f ine 
interpolation operator . 

Three levels of grids are used in the present work. The cor- 
rections to the potential in each grid are computed in the half 
cycle from fine grid to coarse grid, while the results from the 
coarser grid are interpolated to the finer grid in the half cycle 
from coarse grid to fine grid. 

In this preliminary study, no attempt has been made to change 
the method of applying the boundary conditions. They are upoated at 
the eno of the smoothing operations. This turns out to be one of 
the reasons tliat the preliminary computation does not converge. 

Each tinie after the potential on the boundary is updated using the 
boundary condition, high-frequency error components are reintro- 
ducei.''. They are interpolated into the coarse grid as if they were 
the low-frequency components. It was pointed out by Jameson and 
Steinhoff (1981) that the introduction of the relaxing bounoary con- 
ditions causes the oifficulty in applying the multigiid technique in 
the transonic flow problem. A special treatment of the boundary 
condition together with a special formula in defining the resioue 
near the boundary is required for the solution to converge. These 
are some of the focal points for future research. 

There are other sources of difficulty in the cascade problem. 

Ttie mesh system used in the present study has only two upstream 
bounoary cells. As the multigrid cycle goes to the coarser griti, an 
Interior grid point s jch as A in h iqure A becomes a boundary point 
in the coarse grid. Since the freestream conditions are applit.d in 
the boundary cells, the boundary points are governed by uifteretit 
operators than the interior points. This violates tlic assumption 
that the same mesh point is governed by the same type of operator in 
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both the coarse and the fine grids. One possible solution to the 
problem is to modify the mesh system so that at least eight cells 
are located on the upstream boundary. 

In principle, there seems to be no reason why the multigrid 
technique should not work in the computation of transonic flow past 
a cascade. There is certainly some research needed in applying the 
boundary condition correctly in the multigrid computations. 
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6 . Recommendations for the Development of a 

1 hree-Dimensional Code 

From the present investigation, some insight can be obtained in 
the extension of the method to three dimensions. 

Ail mesh generation schemes now in use for cascades are based on 
a conformal transformation of one form or another. If the mesh lines 
wrap aiound the leading edge, as 0-type and C-type meshes do, then the 
orthogonality of the conformal transformation results in a mesh system 
which extends upstream a distance on the order of the pitch distance. 
From the two-dimensional computations, it seems that the flow arouno 
the blade is quite accurately predicted for a pitch-to-chord ratio as 
low as 1. However, if the flow in a long hub ano/or a long shroud 
ahead of the blade is to be included in the computation, some other 
mesh system must be patched upstream. 

To assemble a series of two-dimensional cascade meshes to form a 
three-dimensional mesh system, several additional constraints must be 
satisfied. First of all, the mesh lines leading from the upstream 
boundary to the blade surface on each spanwise two-dimensional mesh 
must be connected so that the three-dimensional upstream boundary can 
be clearly defined. Because of the periodic mesh constraint, the in- 
dex of this line is always the midpoint in the wraparound direction. 
Furthermore, to form a clear-cut trailing edge definition, the trailiny- 
edye indices (ITEl and ITE2) must be maintained for all spanwise 
stations. To accomplish this, a stretching transformation in the 
wraparound direction is requireo. This stretching transformation must 
apply to all pitch-stagger angle combinations in the realistic range. 

For a higfi stagger angle and a large pitch distance, the stretching 
transformation causes most of tfie points to cluster un the upper sui - 
face of the blaue and leaves very few points on the lower surface. 

This is not a completely satisfactory mesh. However, foi the applica- 
tion to turbomachinery, the present mesh generation sctieme is probably 
applicable to all practical ranges of pitch-stagger angle combinations. 
This expectation needs to be qualified by further investigation. 
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Currently, research in the application of the so-called H-type 
mesh to the potential flow calculation is underway (Jou 1960). The 
H-type grid is essentially a streamline-potential line type of grid 
which has a transformation singularity. If the development is 
successful, it may have a profound effect on the development of a 
three-dimensional code around a rotor. 

As to the solution algorithm, the following recommendation is 
given. The artificial density formulation shows a great deal more 
tolerance to the grid irregularity as compared to the artificial 
viscosity. Although the electrostatic mesh generation scheme is 
expected to give a smooth mesh, some irregularity may arise when the 
mesh spacing becomes extremely dense, like near the leading edge. A 
three-dimensional code based on the artificial density formulation 
is preferred over that using artificial viscosity. 
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